{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2004-Particle swarm optimization with particles having quantum behavior\n",
    "\n",
    "## Quantum Delta Potential Well Model of PSO\n",
    "\n",
    "### A. Dynamics of Classical PSO\n",
    "\n",
    "在PSO系统中，有$M$个粒子，每个粒子有自己的位置向量和速度向量（n维空间），表示为：\n",
    "$$ \\vec x_i(t)=\\left(x_{i1}(t),x_{i2}(t),x_{i3}(t),\\cdots, x_{iD}(t)\\right)^\\top $$\n",
    "和\n",
    "$$ \\vec v_i(t)=\\left(v_{i1}(t),v_{i2}(t),v_{i3}(t),\\cdots, v_{iD}(t)\\right)^\\top $$\n",
    "粒子的轨迹公式：\n",
    "$$\n",
    "\\begin{align}\n",
    "\\vec v_i(t+1)&=w\\cdot\\vec v_i(t)+\\varphi_1^\\top(\\vec p_i-\\vec x_i(t))+\\varphi_2^\\top(\\vec p_g-\\vec x_i(t)) \\tag{1a}\n",
    "\\\\\n",
    "\\vec x_i(t+1)&=\\vec x_i+\\vec v_i(t+1) \\tag{1b}\n",
    "\\\\\n",
    "i&=(1,2,\\cdots,M) \\notag\n",
    "\\end{align}\n",
    "$$\n",
    "$\\varphi_1^\\top$和$\\varphi_2^\\top$是随机向量：\n",
    "$$\n",
    "\\begin{align}\n",
    "\\varphi_1^\\top=(\\varphi_{11},\\varphi_{12},\\cdots,\\varphi_{1D}) \\tag{2a}\n",
    "\\\\\n",
    "\\varphi_2^\\top=(\\varphi_{21},\\varphi_{22},\\cdots,\\varphi_{2D}) \\tag{2b}\n",
    "\\\\\n",
    "\\left(\\varphi_{kl}>0,k=1,2\\right) \\notag\n",
    "\\end{align}\n",
    "$$\n",
    "等式$（1a）$中，向量$\\vec p_i=(p_{i1},p_{i2},\\cdots,p_{iD})$是粒子$i$的最好位置，向量$\\vec p_g=(p_{g1},p_{g2},\\cdots,p_{gD})$是粒子群最好位置，参数$w$是惯性权重（PSO中没有），若正确选择$\\varphi_{kl}$，则$\\vec x_i$将会收敛至$\\vec p=(p_1,p_2,\\cdots,p_D)$，其坐标为：\n",
    "$$ p_d=(\\varphi_{1d}p_{id}+\\varphi_{2d}p_{gd})/(\\varphi_{1d}+\\varphi_{2d}) $$\n",
    "为了收敛使用Delta potential well（一种potential field）绑定粒子\n",
    "\n",
    "### B. Quantum Delta Potential Well Model of PSO\n",
    "\n",
    "三维空间中，粒子的波函数$\\Psi(\\vec x,t)$：\n",
    "$$ \\begin{align}|\\Psi|^2dxdydz=Qdxdydz\\tag{4}\\end{align} $$\n",
    "$Qdxdydz$是在$t$时刻测量粒子位置时，发现它位于点$(x,y,z)$周围的体积元中的概率，因此$|\\Psi|^2$为概率密度函数，满足：\n",
    "$$ \\begin{align}\\int_{-\\infty}^{+\\infty}|\\Psi|^2dxdydz=\\int_{-\\infty}^{+\\infty}Qdxdydz\\tag{5}\\end{align} $$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 例子：优化函数\n",
    "\n",
    "假设我们希望最小化一个简单的函数：\n",
    "$$ f(x)=x^2 $$\n",
    "在区间$[−10,10]$内。\n",
    "\n",
    "假设我们初始化5个粒子，每个粒子的位置和速度随机生成：\n",
    "- 粒子位置：$x_1,x_2,x_3,x_4,x_5$\n",
    "- 粒子速度：$v_1,v_2,v_3,v_4,v_5$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在运行过程中，算法会输出每次迭代中最优粒子的适应度值，最终会逐渐收敛到最优解$x=0$对应的函数值为$f(0)=0$。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import random\n",
    "from typing import List\n",
    "\n",
    "\n",
    "class Particle(object):\n",
    "    def __init__(self, bounds) -> None:\n",
    "        self._x = np.zeros(len(bounds))\n",
    "        for idx, (lo, hi) in enumerate(bounds):\n",
    "            self._x[idx] = random.uniform(lo, hi)\n",
    "\n",
    "        self._best = self._x.copy()\n",
    "        self._best_value = np.NaN\n",
    "\n",
    "    def __str__(self) -> str:\n",
    "        return str(self._x)\n",
    "\n",
    "    @property\n",
    "    def best(self):\n",
    "        return self._best\n",
    "\n",
    "    def set_best(self, x) -> None:\n",
    "        for i, v in enumerate(x):\n",
    "            self._best[i] = v\n",
    "\n",
    "    @property\n",
    "    def best_value(self):\n",
    "        return self._best_value\n",
    "\n",
    "    def set_best_value(self, v) -> None:\n",
    "        self._best_value = v\n",
    "\n",
    "    def __getitem__(self, key):\n",
    "        return self._x[key]\n",
    "\n",
    "    def __setitem__(self, key, val) -> None:\n",
    "        self._x[key] = val\n",
    "\n",
    "\n",
    "class Swarm(object):\n",
    "    def __init__(self, size, dim, bounds) -> None:\n",
    "        self._particles = [Particle(bounds) for _ in range(0, size)]\n",
    "        self._dim = dim\n",
    "        self._gbest_value = None\n",
    "        self._gbest = None\n",
    "\n",
    "    def size(self) -> int:\n",
    "        return len(self._particles)\n",
    "\n",
    "    def particles(self) -> List[Particle]:\n",
    "        return self._particles\n",
    "\n",
    "    def mean_best(self):\n",
    "        x = np.zeros(self._dim)\n",
    "        for particle in self._particles:\n",
    "            x += particle.best\n",
    "\n",
    "        x /= self.size()\n",
    "\n",
    "        return x\n",
    "\n",
    "    @property\n",
    "    def gbest(self):\n",
    "        return self._gbest\n",
    "\n",
    "    @property\n",
    "    def gbest_value(self):\n",
    "        return self._gbest_value\n",
    "\n",
    "    def update_gbest(self) -> None:\n",
    "        get_best_value = lambda x: x.best_value\n",
    "        pg = min(self.particles(), key=get_best_value)\n",
    "        if (not self._gbest_value) or self._gbest_value > pg.best_value:\n",
    "            self._gbest = pg.best.copy()\n",
    "            self._gbest_value = pg.best_value\n",
    "\n",
    "\n",
    "class QPSO(Swarm):\n",
    "    def __init__(self, cf, size, dim, bounds, max_iters) -> None:\n",
    "        super(QPSO, self).__init__(size, dim, bounds)\n",
    "        self._cf = cf\n",
    "        self._max_iters = max_iters\n",
    "        self._iters = 0\n",
    "        self.init_eval()\n",
    "\n",
    "    def init_eval(self) -> None:\n",
    "        for particle in self.particles():\n",
    "            best_value = self._cf(particle[:])\n",
    "            particle.set_best_value(best_value)\n",
    "\n",
    "        self.update_gbest()\n",
    "\n",
    "    def update_best(self) -> None:\n",
    "        for particle in self.particles():\n",
    "            best_value = self._cf(particle[:])\n",
    "            if best_value < particle.best_value:\n",
    "                particle.set_best(particle[:])\n",
    "                particle.set_best_value(best_value)\n",
    "\n",
    "        self.update_gbest()\n",
    "\n",
    "    def kernel_update(self, **kwargs) -> None:\n",
    "        pass\n",
    "\n",
    "    def update(self, callback=None, interval=None) -> None:\n",
    "        while self._iters <= self._max_iters:\n",
    "            self.kernel_update()\n",
    "            self.update_best()\n",
    "            if callback and (self._iters % interval == 0):\n",
    "                callback(self)\n",
    "\n",
    "            self._iters += 1\n",
    "\n",
    "    @property\n",
    "    def iters(self):\n",
    "        return self._iters\n",
    "\n",
    "    @property\n",
    "    def max_iters(self):\n",
    "        return self._max_iters\n",
    "\n",
    "\n",
    "class QDPSO(QPSO):\n",
    "    def __init__(self, cf, size, dim, bounds, max_iters, g) -> None:\n",
    "        super(QDPSO, self).__init__(cf, size, dim, bounds, max_iters)\n",
    "        self._g = g\n",
    "\n",
    "    def kernel_update(self, **kwargs) -> None:\n",
    "        for particle in self._particles:\n",
    "            for d in range(0, self._dim):\n",
    "                varphi_1 = random.uniform(0.0, 1.0)\n",
    "                varphi_2 = random.uniform(0.0, 1.0)\n",
    "                u = random.uniform(0.0, 1.0)\n",
    "                p = (varphi_1 * particle.best[d] + varphi_2 * self._gbest[d]) / (\n",
    "                    varphi_1 + varphi_2\n",
    "                )\n",
    "                L = (1 / self._g) * abs(particle[d] - p)\n",
    "                particle[d] = p + (1 if random.random() > 0.5 else -1) * L * np.log(\n",
    "                    1.0 / u\n",
    "                )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Found best position: [ 2.88708574e-05 -4.04175963e-05  6.91056576e-05 -3.45858201e-05\n",
      " -3.29625115e-05  6.42533189e-05 -1.16080762e-04 -2.96035617e-05\n",
      " -9.50431070e-05 -2.22084285e-04]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAGdCAYAAACyzRGfAAAAP3RFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMS5wb3N0MSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8kixA/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAysUlEQVR4nO3de3zU1Z3/8fd3ZjKT60zIPYEA4aKgIiLIzUutUpFa64V1q2sttVYfutiKtLXSVrf9dS2u7bZuW1a3262uWxVrq7ZaxSoqaOUieEXuEC4BEnIhmdwvM+f3R8hIJECGzMx3MvN6Ph7z0Mx8J/Px6IN5e87nnK9ljDECAACIEYfdBQAAgORC+AAAADFF+AAAADFF+AAAADFF+AAAADFF+AAAADFF+AAAADFF+AAAADHlsruATwsGg9q/f7+ysrJkWZbd5QAAgH4wxqixsVElJSVyOI4/txF34WP//v0qLS21uwwAAHAS9u7dq2HDhh33mrgLH1lZWZK6i/d6vTZXAwAA+sPv96u0tDT0PX48cRc+epZavF4v4QMAgEGmPy0TNJwCAICYInwAAICYInwAAICYInwAAICYInwAAICYInwAAICYInwAAICYInwAAICYInwAAICYInwAAICYInwAAICYInwAAICYSurwsa++Vf+1YodaOwJ2lwIAQNKIu7vaxtJP/rpJf/3ogA61dOruOePsLgcAgKSQtDMfxhi9vaNGkvTH9RXqDARtrggAgOSQtOFj28EmHWrplCTVNLVr+aYqmysCACA5JG34WFNe1+vnJ9futakSAACSS/KGj521kqSrzx4qSVq5rVoVh1rsLAkAgKSQlOHDGKO1h2c+/nFKqWaOzpUx0h/WVdhcGQAAiS8pw8eu2hYdbGyX2+nQWaXZunbqcEnS0+v2KhA0NlcHAEBiS8rwsba8e8nlrNJspaY4Nfv0Qg1JT9GBhjat2HrQ5uoAAEhsSRk+1uzsXnKZWpYjSfK4nLr67GGSaDwFACDakjN8HO73mDYqJ/TcdVNLJUmvbT6og/42W+oCACAZJF34qDjUon31rXI5LE0eMST0/JiCLE0ZMUSBoNHT62k8BQAgWpIufPQsuZwx1Kd0d+/T5XsaT5e+s4fGUwAAoiTpwsfaPpZcelw2oVi+tBTtrWvV65tpPAUAIBqSLnysObzTZVrZ0eEjze3UtYd7Px59e1csywIAIGkkVfio8rdpV22LLEuaMvLo8CFJN0wfIYclvbW9RlurGmNcIQAAiS+pwkfPLpfTir3ypqb0ec2wIem65LQiScx+AAAQDUkVPtaGllxyj3vdjeeOlCQ9826FGg7f+RYAAERGUoWPnp0ufTWbHmlqWY7GF3vV1hnU0nf2xKI0AACSRtKEj9qmdm072CRJOucY/R49LMvSjTNHSpIeW7VbXYFgtMsDACBpJE342F3XotwMt04tzFJOhvuE13/xrBLlZLi1r75Vr25i2y0AAJGSNOHj7OFDtO4Hs/TEzdP6dX1qijN05Pojfy+PZmkAACSVpAkfUvdySm6mp9/Xf3n6CDkdltaU12njfn8UKwMAIHkkVfgIV7EvTbPGF0iSXt/C0gsAAJFA+DiBMQWZkroPKAMAAANH+DiBgqxUSVJ1Y7vNlQAAkBgIHyeQn9XdI3KQ8AEAQEQQPk6gIBQ+WHYBACASCB8n0LPsctDfLmOMzdUAADD4ET5OoGfZpb0rqMb2LpurAQBg8CN8nECa26ksj0tS9+wHAAAYGMJHP+R76fsAACBSCB/90NN0ynZbAAAGjvDRD0c2nQIAgIEhfPRDaOajifABAMBAET76IXTQGEesAwAwYISPfijwcsopAACRQvjoh1DPB+EDAIABI3z0QwHLLgAARAzhox96Zj78bV1q6wzYXA0AAIMb4aMfvGkuuV3dQ8VZHwAADAzhox8sy1J+Jk2nAABEAuGjn3p2vFRzxDoAAANC+OgnjlgHACAyCB/9FDpojPABAMCAED76ifu7AAAQGYSPfgqd9UHPBwAAA0L46CeOWAcAIDIIH/3Us+xCwykAAAND+OinnobTmqZ2BYLG5moAABi8wgofixcv1jnnnKOsrCwVFBToyiuv1JYtW3pd09bWpvnz5ys3N1eZmZmaO3euqqqqIlq0HXIz3LIsKWik2mZmPwAAOFlhhY8VK1Zo/vz5Wr16tV555RV1dnbqkksuUXNzc+iaO++8U88//7yefvpprVixQvv379fVV18d8cJjzeV0KDej5wZzhA8AAE6WK5yLly1b1uvnRx99VAUFBVq/fr0uuOACNTQ06H/+53/0xBNP6KKLLpIkPfLIIxo/frxWr16t6dOnR65yGxRkeVTT1K7qJsIHAAAna0A9Hw0NDZKknJwcSdL69evV2dmpWbNmha4ZN26chg8frlWrVvX5O9rb2+X3+3s94lXoiHVmPgAAOGknHT6CwaAWLFigc889V2eccYYkqbKyUm63W9nZ2b2uLSwsVGVlZZ+/Z/HixfL5fKFHaWnpyZYUdZ/cXI6zPgAAOFknHT7mz5+vDRs2aOnSpQMqYNGiRWpoaAg99u7dO6DfF02c9QEAwMCF1fPR4/bbb9cLL7yglStXatiwYaHni4qK1NHRofr6+l6zH1VVVSoqKurzd3k8Hnk8npMpI+Y4Yh0AgIELa+bDGKPbb79dzz77rF577TWVlZX1en3y5MlKSUnR8uXLQ89t2bJFe/bs0YwZMyJTsY1Cd7al4RQAgJMW1szH/Pnz9cQTT+jPf/6zsrKyQn0cPp9PaWlp8vl8uummm7Rw4ULl5OTI6/XqG9/4hmbMmDHod7pIR97ZNvyej65AUB9U1Ous0iFyOqxIlwYAwKAR1szHQw89pIaGBl144YUqLi4OPZ566qnQNb/4xS/0hS98QXPnztUFF1ygoqIiPfPMMxEv3A5HLrsYE94ppz98/mPNfWiVnly7JxqlAQAwaIQ189GfL9zU1FQtWbJES5YsOemi4lVPw2l7V1D+ti750lL69b7ymmY9uba7kXbVzlp9efqIqNUIAEC8494uYUhNcSortTuvVYex9PKLV7aG7gezaX/8nmMCAEAsED7CVJAV3nbbzZV+Pf/h/tDP5bXNam7vikptAAAMBoSPMPU0nVb3M3z8+9+2yhjpsgnFKsjyyBhpc2VjNEsEACCuET7CFM5ZH+/vrdcrG6vksKQ7P3eKTivxSpI2HWDpBQCQvAgfYSoIY7vtv/9tiyTpqknDNKYgU6cVd4ePjYQPAEASI3yEKXRzuRMsu6zaUas3t9UoxWlpwayxkhSa+dhI0ykAIIkRPsKU34+GU2OMfnZ41uPac4arNCddkjT+8MzH5kp/aPcLAADJhvARplDPx3HCx/++vUvrdx+Sx+XQ7ReNCT0/MjdDaSlOtXUGVV7THPVaAQCIR4SPMPX0fOw71KrKhqP7Pjbu9+snL26WJC2aM06F3tTQa06HpXHFWZJoOgUAJC/CR5hG5mVoVF6GWjsD+qffru7V+9HS0aXbn3xXHYGgZo0v0LyZI496P02nAIBkR/gIU4rTocdumqqh2WnaWd2sL/92jeqaOyRJP/rLRu2sblah16MH/mGiLOvoG8jRdAoASHaEj5MwbEi6Hv/6NBV6PdpS1agv/3aNHl+zW0+t2yvLkn7xpbOUk+Hu873jmfkAACQ5wsdJGpmXoce/Pl15mW5tPODX95/dIEm6/bNjNHN03jHfN64oS5bVvVW3v6ekAgCQSAgfAzCmIFOPf326hqR339327OHZuuPiscd9T7rbpbK8DEk0nQIAkhPhY4BOLcrS07fO0DcvGqOHb5gsl/PEQ0rTKQAgmRE+ImBMQZYWXnJq6AyQE6HpFACQzAgfNmDmAwCQzAgfNugJHzurm9TWGbC5GgAAYovwYYP8LI/yMt0KGmlLZaPd5QAAEFOEDxtYlsV5HwCApEX4sAlNpwCAZEX4sAlNpwCAZEX4sMnph2c+Nh3wyxhjczUAAMQO4cMmpTnpkqSWjoD8rV02VwMAQOwQPmzicTmVleqSJFU3cY8XAEDyIHzYKD/TI0mqIXwAAJII4cNGuZluSVJtU4fNlQAAEDuEDxvlMfMBAEhChA8bfTLzQfgAACQPwoeNemY+qll2AQAkEcKHjXrCBzMfAIBkQviwUd7hZRd6PgAAyYTwYaNPGk5ZdgEAJA/Ch41yWXYBACQhwoeNepZdmjsCau0I2FwNAACxQfiwUabHJber+18BfR8AgGRB+LCRZVkcsQ4ASDqED5vlccQ6ACDJED5slsvMBwAgyRA+bMZZHwCAZEP4sFkuZ30AAJIM4cNm3NkWAJBsCB82o+EUAJBsCB82Y+YDAJBsCB82C93ZtpmZDwBAciB82Cz38LLLoZYOdQWCNlcDAED0ET5sNiTdLYclGSPVMfsBAEgChA+bOR2WcjJ6zvogfAAAEh/hIw7QdAoASCaEjzjwSdMp4QMAkPgIH3Ggp+m0ppFlFwBA4iN8xIHQsgszHwCAJED4iAPMfAAAkgnhIw7QcAoASCaEjzgQur8Lyy4AgCRA+IgDoZkPll0AAEmA8BEHjtxqa4yxuRoAAKKL8BEHek447QwY+Vu7bK4GAIDoInzEgdQUp7JSXZLYbgsASHyEjzjxSd8H4QMAkNjCDh8rV67U5ZdfrpKSElmWpeeee67X61/96ldlWVavx6WXXhqpehNWz44Xbi4HAEh0YYeP5uZmTZw4UUuWLDnmNZdeeqkOHDgQejz55JMDKjIZcH8XAECycIX7hjlz5mjOnDnHvcbj8aioqOiki0pGn5xySvgAACS2qPR8vPHGGyooKNCpp56q2267TbW1tce8tr29XX6/v9cjGX1yfxeWXQAAiS3i4ePSSy/VY489puXLl+vf/u3ftGLFCs2ZM0eBQKDP6xcvXiyfzxd6lJaWRrqkQSGXhlMAQJIIe9nlRK699trQ30+YMEFnnnmmRo8erTfeeEMXX3zxUdcvWrRICxcuDP3s9/uTMoDkh45YZ+YDAJDYor7VdtSoUcrLy9P27dv7fN3j8cjr9fZ6JKNcbi4HAEgSUQ8fFRUVqq2tVXFxcbQ/alDjnA8AQLIIe9mlqamp1yxGeXm53n//feXk5CgnJ0c/+tGPNHfuXBUVFWnHjh266667NGbMGM2ePTuihSeannM+mjsCau0IKM3ttLkiAACiI+yZj3Xr1mnSpEmaNGmSJGnhwoWaNGmS7r33XjmdTn344Yf64he/qFNOOUU33XSTJk+erDfffFMejyfixSeSTI9Lblf3vw6WXgAAiSzsmY8LL7zwuHdeffnllwdUULKyLEv5mR7tq29VbXOHSnPS7S4JAICo4N4ucYSDxgAAyYDwEUc4Yh0AkAwIH3Gkp+m0soHwAQBIXISPOHJ6iU+S9NeP9h+3rwYAgMGM8BFHrjp7qNLdTm2tatKqHce+Hw4AAIMZ4SOOeFNTNPfsYZKkR9/eZW8xAABECeEjzsybOUKS9OqmKu2ta7G5GgAAIo/wEWfGFGTp/LF5Chrp96t3210OAAARR/iIQ/NmjJQkLX1nr1o7AvYWAwBAhBE+4tBnxxVoeE66Glo79dz7++wuBwCAiCJ8xCGnw9JXZnT3fvzv27vYdgsASCiEjzh1zZRSpaU4tbmyUat31vV5TSBo9NJHB3T9b1frjqXvKRAkpAAA4l/YN5ZDbPjSUnT12UP1+Jo9+sUrW9XQWqaReekanpOuQNDoD+sq9Ojb5dpb1xp6z0XjCnTFWUNtrBoAgBMjfMSxr84cqcfX7NHaXXVau+uT2Y8Up6XOQPcsx5D0FJ1W4tXft9fqF69s1ecnFCvFyYQWACB+ET7i2NjCLP37NRO1Ymu1dtc2a3ddi+pbOtUZMBqdn6Gbzhulq88eqkDQ6IIHXteu2hY9826FvnTOcLtLBwDgmCwTZ92Mfr9fPp9PDQ0N8nq9dpcTdxpaOlXX0qEROelyOKzQ8799c6f+9a+bNDQ7Ta99+zPyuJw2VgkASDbhfH8zPz/I+NJTVJaX0St4SNKXp49QodejffWtWrp2r03VAQBwYoSPBJGa4tQ3LhorSfr169s5nAwAELcIHwnkH6eUatiQNFU3tuuxVbvsLgcAgD4RPhKI2+XQHRd3z348tGKHNlf6tbu2WfvqW3XQ36b2LmZDAAD2o+E0wXQFgrrkwZXaWd181GupKQ7NHJ2nz56ar8+OK9CwId1nhpTXNGtzpV9bKhs1Oj9TV07irBAAQHjC+f4mfCSglVurdfefPlRzR0BdgaA6g0ZdgaA+fQDq0Ow01TS1q70rGHrOsqS135ul/CxPjKsGAAxm4Xx/c85HArrglHy9vejiXs8ZY7SlqlGvbT6oNzZXa/2eQ9pX3306alqKU6cWZWnHwSY1tndpR3UT4QMAEDWEjyRhWZbGFXk1rsirf75wjBpaOrVhf4NKstM0PCddToeleb9bqxVbq1Ve06zpo3LtLhkAkKBoOE1SvvQUnTsmT2V5GXIePjOkLC9DklRec3S/CAAAkUL4QMio/O7w0VezKgAAkUL4QMgnMx9NNlcCAEhkhA+E9ISPPXUtCnx6awwAABFC+EBIiS9NbpdDnQGjfYda7S4HAJCgCB8IcTgsleUe7vtg6QUAECWED/QyMi9dEjteAADRQ/hAL2V5mZIIHwCA6CF8oJdRnPUBAIgywgd6KeOsDwBAlBE+0EvPdtv9Da1q6wzYXA0AIBERPtBLboZbWakuGdN93gcAAJFG+EAvlmWF+j5YegEARAPhA0fhBnMAgGgifOAon2y35aAxAEDkET5wlJ4dL8x8AACigfCBo/QcsU74AABEA+EDR+k5Yr2mqUMNrZ02VwMASDSEDxwlKzVF+VkeSdIuZj8AABFG+ECfena87KolfAAAIovwgT5x1gcAIFoIH+gTZ30AAKKF8IE+ET4AANFC+ECfRh1x1ocxxuZqAACJhPCBPpXmpMthSU3tXapuare7HABAAiF8oE8el1NDh6RJksppOgUARBDhA8fUc4+XbQe5xwsAIHJcdheA+DUqL0Mrt1brB89t0OIXN6nQm6r8LI/OG5Onb1w81u7yAACDFDMfOKbLzixWboZbktTcEdDOmmatKa/Tv7+yVXtqW2yuDgAwWDHzgWM6Z2SO1t/zOTW3d+lgY7uq/G363jMfaWdNszZX+jU8N93uEgEAgxAzHzihDI9LZXkZmj4qV2eVZkuStlY12lsUAGDQInwgLKcUZUmStlTRhAoAODmED4Tl1MLu8LG1kpkPAMDJIXwgLD0zHzuqm9TRFbS5GgDAYET4QFhKfKnK9LjUFTTaVcvhYwCA8BE+EBbLsnRKYffhY1tYegEAnISww8fKlSt1+eWXq6SkRJZl6bnnnuv1ujFG9957r4qLi5WWlqZZs2Zp27ZtkaoXceDUw0sv7HgBAJyMsMNHc3OzJk6cqCVLlvT5+gMPPKBf/vKXevjhh7VmzRplZGRo9uzZamtrG3CxiA+nHG46ZeYDAHAywj5kbM6cOZozZ06frxlj9OCDD+oHP/iBrrjiCknSY489psLCQj333HO69tprB1Yt4kJoxwszHwCAkxDRno/y8nJVVlZq1qxZoed8Pp+mTZumVatW9fme9vZ2+f3+Xg/Et54dL7vrWtTaEbC5GgDAYBPR8FFZWSlJKiws7PV8YWFh6LVPW7x4sXw+X+hRWloayZIQBXmZHuVmuGWMtJ073gIAwmT7bpdFixapoaEh9Ni7d6/dJaEfQn0fLL0AAMIU0fBRVFQkSaqqqur1fFVVVei1T/N4PPJ6vb0eiH/seAEAnKyIho+ysjIVFRVp+fLloef8fr/WrFmjGTNmRPKjYDN2vAAATlbYu12ampq0ffv20M/l5eV6//33lZOTo+HDh2vBggX613/9V40dO1ZlZWW65557VFJSoiuvvDKSdcNmpxZ1HzTGzAcAIFxhh49169bps5/9bOjnhQsXSpLmzZunRx99VHfddZeam5t1yy23qL6+Xuedd56WLVum1NTUyFUN2409PPNxoKFNDa2d8qWl2FwRAGCwsIwxxu4ijuT3++Xz+dTQ0ED/R5ybuXi59je06Y+3ztCUkTl2lwMAsFE439+273bB4NVz3gc7XgAA4SB84KSFTjql6RQAEAbCB04aZ30AAE4G4QMnreesjy2VjYqz1iEAQBwjfOCkjSnIlMOSDrV0qrqp3e5yAACDBOEDJy01xamRuRmSpK2V3OMFANA/hA8MCH0fAIBwET4wID19Hy9+dEDBIH0fAIATI3xgQK6ZMkzpbqfW7z6k36/ZbXc5AIBBgPCBARk2JF3fvXScJOnfXtqsikMtNlcEAIh3hA8M2A3TR+ickUPU3BHQomc+YtstAOC4CB8YMIfD0v1zz5Tb5dCb22r0x/UVdpcEAIhjhA9ExOj8TN056xRJ0o9f2KiD/jabKwIAxCvCByLm5vPLNGGoT/62Lt31pw91sJEAAgA4GuEDEeNyOvTAP5wpl8PSG1uqNXPxa7r1/9Zr5dZqtuECAEIIH4io8cVe/XbeFE0eMURdQaNlH1fqK79bq8/87HV9VNFgd3kAgDhA+EDEXXhqgf5020wtW3C+vjpzpLJSXdpb16qfvLjJ7tIAAHGA8IGoGVfk1Q+/eLr++o3zJUlrd9XpUHOHzVUBAOxG+EDUDc9N1/hirwJBo1c3VdldDgDAZoQPxMTs0wslSS9/TPgAgGRH+EBMXHJakSTpzW3VaunosrkaAICdCB+IifHFWSrNSVN7V1Art1bbXQ4AwEaED8SEZVmafXj2g6UXAEhuhA/EzOwzusPH8k1V6gwEba4GAGAXwgdi5uzhQ5Sb4Za/rUurd9baXQ4AwCaED8SM02Hpc6d173r5G0svAJC0CB+Iqdmndy+9/G1jJfd7AYAkRfhATM0YnasMt1NV/nZ9UFFvdzkAABsQPhBTqSlOXTiuQFLvXS9dgaCq/G3MhgBAEnDZXQCSz+zTi/TXDw/oT+9WaEd1k3ZWN2lPXYs6A0Z5mR5dNC5fs8YX6ryxeUp3858oACQayxgTV/+r6ff75fP51NDQIK/Xa3c5iILGtk5N/vGr6jjBdluPy6HZpxfpvqvOUFZqSoyqAwCcjHC+v/nfSsRcVmqKfnndWVpTXqeRuRkalZ+hsrwM5WV6tG7XIb26qUqvbqpSxaFW/eWD/dp7qEWP3jhVvjQCCAAkAmY+EJeMMVq3+5Bufmyd6ls6deYwnx772lRlp7vtLg0A0Idwvr9pOEVcsixL54zM0RNfn66cDLc+rGjQP/33GtU1d9hdGgBggAgfiGunlXj15M3TlZfp1sYDfl33m9WqONRid1kAgAEgfCDunVqUpaW3zFBBlkdbqhp1wQOv6+v/+46Wb6pSgK25ADDo0POBQaO8plnfe+YjrTrivjDFvlTdeO5I3Xz+KFmWZWN1AJDc2O2ChFSWl6Enb5muHdVNenLNHv3x3QodaGjTT17cLIdl6evnj7K7RABAP7DsgkFndH6mfvCF07R60cX61udOkSTd/9JmvbfnkM2VAQD6g/CBQSs1xanbLxqjyyYUqytodPsT76mhpdPusgAAJ0D4wKBmWZYWz52gEbnp2lffqm//8QPFWRsTAOBTCB8Y9LypKVryT2fL7XTolY1V+t3fd9ldEgDgOAgfSAhnDPXpB18YL0m6/6VNenjFDr2ysUpbKhvV0tFlc3UAgCOx2wUJ44bpI7R6Z61e/KhS97+0uddrZw7z6bfzpqggK9Wm6gAAPZj5QMKwLEs/u2aivn3JKfr8hCKdMdQbuhndhxUN+tYfPlCQQ8kAwHbMfCChpLtduv2isb2e23TAr6v+8+96c1uNfvPmTt36mdE2VQcAkJj5QBIYX+zVDy8/XZL0s5e3cB4IANiM8IGk8KVzSnXZmd3ngXxz6Xvyt3EeCADYhfCBpGBZln5y1QQNG5KmvXWt+v6zGzgPBABsQs8HkoYvLUX/ce0k/eN/rdLzH+zX1spGdQaDausIqLUzoEJvqv7vpmnKz/LYXSoAJDRmPpBUJo8Yom9d0n0/mC1VjdpZ3az9DW061NKpzZWN+uFfPra5QgBIfMx8IOnc9pnRmlQ6RO1dAaWlOJXmdqq2uUNf/991+utHB3TFx5W65PQiu8sEgIRF+EDSsSxLM0bnHvX8LReM0kNv7NA9f96g6aNz5U1NsaE6AEh8LLsAh91x8ViV5WWoyt+uxS9uPvEbAAAnhfABHJaa4tT9V0+QJD25do9W7ai1uSIASEyED+AI00bl6p+mDZckLXrmQ7V1BmyuCAASDz0fwKfcPWeclm+q0q7aFp3/wOsq8aUqP8uj/CyPyvIydNWkYWzHBYABsEycnbTk9/vl8/nU0NAgr9drdzlIUq9vOajbfr9ebZ3Bo15zOx267MxizZs5UmeVZse+OACIQ+F8fxM+gGOoa+7QnroWVTe2q7qxXQcb27Ria7Xe21MfumZiabbu/cJ4TR6RY1+hABAHCB9AFH1YUa9H396lFz44oI5AUFmpLv15/rkalZ9pd2kAYJtwvr8j3nD6wx/+UJZl9XqMGzcu0h8D2ObMYdn6+T+epVWLLtI5I4eosa1LNz+2To3crA4A+iUqu11OP/10HThwIPR46623ovExgK1yMz36z+snq8ibqh3VzbrzqQ8UDMbVRCIAxKWohA+Xy6WioqLQIy8vLxofA9guP8uj/7phstwuh17dVKUHl2+zuyQAiHtRCR/btm1TSUmJRo0apeuvv1579uw55rXt7e3y+/29HsBgMrE0W4uv6j6c7JfLt2nZhkqbKwKA+Bbx8DFt2jQ9+uijWrZsmR566CGVl5fr/PPPV2NjY5/XL168WD6fL/QoLS2NdElA1M2dPExfO7dMkvStP7yvg41tNlcEAPEr6rtd6uvrNWLECP385z/XTTfddNTr7e3tam9vD/3s9/tVWlrKbhcMOl2BoC55cKV2Vjfrv26YrNncGRdAEglnt0vUTzjNzs7WKaecou3bt/f5usfjkcfDaZEY/FxOh8YXe7Wzull761rsLgcA4lbU7+3S1NSkHTt2qLi4ONofBdhu2JA0SVLFoVabKwGA+BXx8PHtb39bK1as0K5du/T222/rqquuktPp1HXXXRfpjwLiTumQdElSxSFmPgDgWCK+7FJRUaHrrrtOtbW1ys/P13nnnafVq1crPz8/0h8FxJ3SnO7wsbeOmQ8AOJaIh4+lS5dG+lcCg0bp4WWXvYdaZIyRZVk2VwQA8SfqPR9AMhl6OHy0dARU19xhczUAEJ8IH0AEeVxOFXq7d2/tpekUAPpE+AAijKZTADg+wgcQYTSdAsDxET6ACDuy6RQAcDTCBxBhw4b0zHwQPgCgL4QPIMKG5XDKKQAcD+EDiLCehtN9h1oVDEb1vo0AMCgRPoAIK/alyumw1BEI6mBj+4nfAABJhvABRJjL6VBJdqokmk4BoC+EDyAKhmXTdAoAx0L4AKKg9HDTKWd9AMDRCB9AFHDKKQAcG+EDiILQKaeEDwA4CuEDiAKWXQDg2AgfQBT0nHJ6oKFVnYGgzdUAQHwhfABRkJ/pkdvlUNBIB+rb7C4HAOIK4QOIAofD0rAhPces0/cBAEcifABR0rPjhaZTAOiN8AFECU2nANA3wgcQJcOY+QCAPhE+gCgJLbtwxDoA9EL4AKIktOxyiGUXADgS4QOIkp6Zj+rGdrV1BmyuBgDiB+EDiJLs9BRlelySpApmPwAghPABRIllfXLWB02nAPAJwgcQRT07XipoOgWAEMIHEEU0nQLA0QgfQBT1NJ3uqWXmAwB6ED6AKBpXnCVJWrG1WtWN7TZXAwDxgfABRNGMUbmaWJqt1s6Afv3aNrvLAYC4QPgAosiyLH139qmSpCfW7uG0UwAQ4QOIuplj8nTemDx1Box+8epWu8sBANsRPoAY+M7h2Y9n39unLZWNNlcDAPYifAAxMLE0W3POKJIx0s/+tsXucgDAVoQPIEa+dckpcljSKxur9O6eQ3aXAwC2IXwAMTKmIEtzzx4mSXpg2WZ1BYI2VwQA9iB8ADG04HOnyO10aPXOOk2571V95+kP9NrmKrV3cddbAMnDMsYYu4s4kt/vl8/nU0NDg7xer93lABH3h3V7df9Lm1XX3BF6Lsvj0vXTR2jBrLFKTXHaWB0AnJxwvr8JH4ANugJBrd1Vp2UbKrVsQ6UOHj79tCwvQ/dfPUHTRuXaXCEAhIfwAQwiwaDR3zZW6V/+skFV/u4QcsP0EfrunHHK9Lhsrg4A+iec7296PgCbORyWLj2jSH+78zO6bmqpJOn/Vu/WpQ+uVE0T94MBkHgIH0Cc8KWlaPHVZ+qJr09TiS9VFYda9dQ7e+0uCwAijvABxJmZY/K04HOnSJKeebdCcbYyCgADRvgA4tCcM4rkcTm0o7pZG/b57S4HACKK8AHEoazUFH3utEJJ3feDAYBEQvgA4tTVZw+VJP3lg/2chgogoRA+gDh1/th85Wa4VdPUrje319hdDgBEDOEDiFMpTocun1giSXqOpRcACYTwAcSxqyZ1L728/HGlmtq7bK4GACKD8AHEsTOH+TQqL0NtnUEt21BpdzkAEBGEDyCOWZYVmv1g6QVAoiB8AHHuysPh4+87alTZ0GZzNQAwcNy1CohzpTnpmjoyR2t31emnL2/RhafmKyvVJW9aioaku1XsS1VqitPuMgGg3wgfwCBw5aShWrurTn96t0J/erfiqNfzMt0amp2mkuw0pbm7g4glS5bV/boxkpGRDp/UPjIvQ1PLcnRWaTbBBUDMET6AQWDu5KEqr2nS7toW+ds65W/tUmN7p2qbOtTSEVBNU4dqmjr0QUVDWL/X7XRoYqlPE4Zmy5PikMOSHJYlR09qOYLLYSkn0628TI/yMj3Kz/QoJ9OtDLdTVh/XA8CxWCbO7lrl9/vl8/nU0NAgr9drdzlAXDPGqKG1U/vqW7XvUKv217eqIxA8PNPRc41kWVJPPAgYo437/VpbXqeDje0DrsHtcig3w62cjO5gMjwnXSNy0zUyN0MjctNVmpPO7AqQBML5/mbmAxjELMtSdrpb2elunV7iC+u9xhjtrm3R2vI6ba9uUlfAKGi6H4Gg0ZGTGZYsdQaCqmnqUHVTu2oa21XT1K72rqA6uoI60NCmA8dphs3L9GjYkDQNHZKm0XkZumZKqUpz0k/2HxvAIMfMB4CTYoxRS0dAdc0d3Y+WDh30t2lPXYt21bZod22zdte0qLGPw9GcDkuXn1ms2y4co1OLsmyoHkCkhfP9HbXwsWTJEv30pz9VZWWlJk6cqF/96leaOnXqCd9H+AAShzFG/tYu7T3UoopDrao41KIVW6v15rZP7lVz8bgCXTOlVGeVZqvQ66F/BBikbA8fTz31lL7yla/o4Ycf1rRp0/Tggw/q6aef1pYtW1RQUHDc9xI+gMS3YV+DHnpjh17ccEBH/gmUn+XRmUN9Gl/slTfNpXS3Sxkep9LdLvnSUkK9JdnpbjkdhBQgntgePqZNm6ZzzjlHv/71ryVJwWBQpaWl+sY3vqG77777uO8lfADJY2d1kx59e5fWltdp28EmBYL9++PIYUm+tBRlpaYo0+NSZqpLmR6X3E6HHI5Pthk7rD7++qnf5XJaKshKVUl2qkqy01TsS5M3zSWHZclpWXI4LLkclpwOSylOB6EHOAZbG047Ojq0fv16LVq0KPScw+HQrFmztGrVqqOub29vV3v7Jx33fr8/0iUBiFOj8jP1/644Q5LU2hHQxgN+fVRRr+3VTWpuD6ilo0stHQE1tXepoaVTtc0damjtVNBIh1o6dailM+Y1W5aU4ugOOUe9dlS0ObEUp6Ws1BR501KUdThE9bXV+Xj1SN2BzOV0hIKSy2EdDl1WaAv1kbuejv87CVixYOcw52V6NP+zY2z7/IiHj5qaGgUCARUWFvZ6vrCwUJs3bz7q+sWLF+tHP/pRpMsAMMikuZ2aPGKIJo8YctzrOgNBHWrpUH1LpxrbutTU3qWmti41tXeqI2AkYxQ03f0mgcN/NUbdu3j6mOjt6Aqqyt+m/fVtOtDQqgP1bWrq6NKx5oSNkToCQSkQiX9qqbVT8rd1aV99a2R+IdAPo/IzEit8hGvRokVauHBh6Ge/36/S0lIbKwIQz1KcDhVkpaogKzWqn2MOh5hAsHvrcVcwqM6AUVcgqM6gUbCfS0Qn0t4VVGNbd5BqPByi+rsYbqTQtQHTXVNX0ChwuNaef4bgEYEsnny6HKP4qi+RDUl32/r5EQ8feXl5cjqdqqqq6vV8VVWVioqKjrre4/HI4/FEugwAGBDLsuS0dESPBwelAZES8bvaut1uTZ48WcuXLw89FwwGtXz5cs2YMSPSHwcAAAaZqCy7LFy4UPPmzdOUKVM0depUPfjgg2pubtaNN94YjY8DAACDSFTCx5e+9CVVV1fr3nvvVWVlpc466ywtW7bsqCZUAACQfDheHQAADFg4398R7/kAAAA4HsIHAACIKcIHAACIKcIHAACIKcIHAACIKcIHAACIKcIHAACIKcIHAACIKcIHAACIqagcrz4QPQeu+v1+mysBAAD91fO93Z+D0+MufDQ2NkqSSktLba4EAACEq7GxUT6f77jXxN29XYLBoPbv36+srCxZlhXR3+33+1VaWqq9e/dy35goYpxjg3GODcY5dhjr2IjWOBtj1NjYqJKSEjkcx+/qiLuZD4fDoWHDhkX1M7xeL/9hxwDjHBuMc2wwzrHDWMdGNMb5RDMePWg4BQAAMUX4AAAAMZVU4cPj8ehf/uVf5PF47C4loTHOscE4xwbjHDuMdWzEwzjHXcMpAABIbEk18wEAAOxH+AAAADFF+AAAADFF+AAAADGVNOFjyZIlGjlypFJTUzVt2jStXbvW7pIGtcWLF+ucc85RVlaWCgoKdOWVV2rLli29rmlra9P8+fOVm5urzMxMzZ07V1VVVTZVnBjuv/9+WZalBQsWhJ5jnCNn3759+vKXv6zc3FylpaVpwoQJWrduXeh1Y4zuvfdeFRcXKy0tTbNmzdK2bdtsrHjwCQQCuueee1RWVqa0tDSNHj1aP/7xj3vdD4RxDt/KlSt1+eWXq6SkRJZl6bnnnuv1en/GtK6uTtdff728Xq+ys7N10003qampKToFmySwdOlS43a7ze9+9zvz8ccfm5tvvtlkZ2ebqqoqu0sbtGbPnm0eeeQRs2HDBvP++++bz3/+82b48OGmqakpdM2tt95qSktLzfLly826devM9OnTzcyZM22senBbu3atGTlypDnzzDPNHXfcEXqecY6Muro6M2LECPPVr37VrFmzxuzcudO8/PLLZvv27aFr7r//fuPz+cxzzz1nPvjgA/PFL37RlJWVmdbWVhsrH1zuu+8+k5uba1544QVTXl5unn76aZOZmWn+4z/+I3QN4xy+F1980Xz/+983zzzzjJFknn322V6v92dML730UjNx4kSzevVq8+abb5oxY8aY6667Lir1JkX4mDp1qpk/f37o50AgYEpKSszixYttrCqxHDx40EgyK1asMMYYU19fb1JSUszTTz8dumbTpk1Gklm1apVdZQ5ajY2NZuzYseaVV14xn/nMZ0Lhg3GOnO9+97vmvPPOO+brwWDQFBUVmZ/+9Keh5+rr643H4zFPPvlkLEpMCJdddpn52te+1uu5q6++2lx//fXGGMY5Ej4dPvozphs3bjSSzDvvvBO65qWXXjKWZZl9+/ZFvMaEX3bp6OjQ+vXrNWvWrNBzDodDs2bN0qpVq2ysLLE0NDRIknJyciRJ69evV2dnZ69xHzdunIYPH864n4T58+frsssu6zWeEuMcSX/5y180ZcoUXXPNNSooKNCkSZP03//936HXy8vLVVlZ2WusfT6fpk2bxliHYebMmVq+fLm2bt0qSfrggw/01ltvac6cOZIY52joz5iuWrVK2dnZmjJlSuiaWbNmyeFwaM2aNRGvKe5uLBdpNTU1CgQCKiws7PV8YWGhNm/ebFNViSUYDGrBggU699xzdcYZZ0iSKisr5Xa7lZ2d3evawsJCVVZW2lDl4LV06VK9++67euedd456jXGOnJ07d+qhhx7SwoUL9b3vfU/vvPOOvvnNb8rtdmvevHmh8ezrzxLGuv/uvvtu+f1+jRs3Tk6nU4FAQPfdd5+uv/56SWKco6A/Y1pZWamCgoJer7tcLuXk5ERl3BM+fCD65s+frw0bNuitt96yu5SEs3fvXt1xxx165ZVXlJqaanc5CS0YDGrKlCn6yU9+IkmaNGmSNmzYoIcffljz5s2zubrE8Yc//EGPP/64nnjiCZ1++ul6//33tWDBApWUlDDOSSThl13y8vLkdDqP6v6vqqpSUVGRTVUljttvv10vvPCCXn/9dQ0bNiz0fFFRkTo6OlRfX9/resY9POvXr9fBgwd19tlny+VyyeVyacWKFfrlL38pl8ulwsJCxjlCiouLddppp/V6bvz48dqzZ48khcaTP0sG5jvf+Y7uvvtuXXvttZowYYJuuOEG3XnnnVq8eLEkxjka+jOmRUVFOnjwYK/Xu7q6VFdXF5VxT/jw4Xa7NXnyZC1fvjz0XDAY1PLlyzVjxgwbKxvcjDG6/fbb9eyzz+q1115TWVlZr9cnT56slJSUXuO+ZcsW7dmzh3EPw8UXX6yPPvpI77//fugxZcoUXX/99aG/Z5wj49xzzz1qu/jWrVs1YsQISVJZWZmKiop6jbXf79eaNWsY6zC0tLTI4ej91eN0OhUMBiUxztHQnzGdMWOG6uvrtX79+tA1r732moLBoKZNmxb5oiLewhqHli5dajwej3n00UfNxo0bzS233GKys7NNZWWl3aUNWrfddpvx+XzmjTfeMAcOHAg9WlpaQtfceuutZvjw4ea1114z69atMzNmzDAzZsywserEcORuF2MY50hZu3atcblc5r777jPbtm0zjz/+uElPTze///3vQ9fcf//9Jjs72/z5z382H374obniiivYAhqmefPmmaFDh4a22j7zzDMmLy/P3HXXXaFrGOfwNTY2mvfee8+89957RpL5+c9/bt577z2ze/duY0z/xvTSSy81kyZNMmvWrDFvvfWWGTt2LFttB+pXv/qVGT58uHG73Wbq1Klm9erVdpc0qEnq8/HII4+ErmltbTX//M//bIYMGWLS09PNVVddZQ4cOGBf0Qni0+GDcY6c559/3pxxxhnG4/GYcePGmd/85je9Xg8Gg+aee+4xhYWFxuPxmIsvvths2bLFpmoHJ7/fb+644w4zfPhwk5qaakaNGmW+//3vm/b29tA1jHP4Xn/99T7/TJ43b54xpn9jWltba6677jqTmZlpvF6vufHGG01jY2NU6rWMOeJYOQAAgChL+J4PAAAQXwgfAAAgpggfAAAgpggfAAAgpggfAAAgpggfAAAgpggfAAAgpggfAAAgpggfAAAgpggfAAAgpggfAAAgpggfAAAgpv4/eJz3Ed160PYAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "\n",
    "def sphere(args):\n",
    "    f = sum([np.power(x, 2.0) for x in args])\n",
    "    return f\n",
    "\n",
    "\n",
    "y = []\n",
    "\n",
    "\n",
    "def log(s):\n",
    "    best_value = [particle.best_value for particle in s.particles()]\n",
    "    best_value_avg = np.mean(best_value)\n",
    "    best_value_std = np.std(best_value)\n",
    "    # print(\n",
    "    #     \"{0: >5}  {1: >9}  {2: >9}  {3: >9}\".format(\n",
    "    #         \"Iters.\", \"Best\", \"Best(Mean)\", \"Best(STD)\"\n",
    "    #     )\n",
    "    # )\n",
    "    # print(\n",
    "    #     \"{0: >5}  {1: >9.3E}  {2: >9.3E}  {3: >9.3E}\".format(\n",
    "    #         s.iters, s.gbest_value, best_value_avg, best_value_std\n",
    "    #     )\n",
    "    # )\n",
    "    y.append(best_value_std)\n",
    "\n",
    "\n",
    "NParticle = 40\n",
    "MaxIters = 100\n",
    "NDim = 10\n",
    "bounds = [(-2.56, 5.12) for i in range(0, NDim)]\n",
    "g = 0.96\n",
    "s = QDPSO(sphere, NParticle, NDim, bounds, MaxIters, g)\n",
    "s.update(callback=log, interval=1)\n",
    "print(\"Found best position: {0}\".format(s.gbest))\n",
    "\n",
    "sns.lineplot(y=y, x=[_ for _ in range(len(y))])\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
